function joint_pos= inverseKinematics(x, y, l1, l2)
    % Compute the inverse kinematics of a 2-R robot
    % Inputs:
    %   x, y: Desired end effector position
    %   l1, l2: Lengths of the robot links
    % Outputs:
    %   q1, q2: Joint angles that achieve the desired position

    % Compute q2 using the law of cosines
    c2 = (x^2 + y^2 - l1^2 - l2^2) / (2*l1*l2);
    s2 = sqrt(1 - c2^2);
    q2 = atan2(s2, c2);

    % Compute q1 using trigonometry
    s1 = ((l1 + l2*c2)*y - l2*s2*x) / (x^2 + y^2);
    c1 = ((l1 + l2*c2)*x + l2*s2*y) / (x^2 + y^2);
    q1 = atan2(s1, c1);

    % Check joint limits
    if q1 < -pi/2
        q1 = q1 + pi;
        q2 = q2 - pi;
    elseif q1 > pi/2
        q1 = q1 - pi;
        q2 = q2 - pi;
    end

    % Normalize joint angles to [-pi, pi]
    q1 = mod(q1 + pi, 2*pi) - pi;
    q2 = mod(q2 + pi, 2*pi) - pi;
    joint_pos = [q1;q2];
end
